function dPds = gradPs(alpha,H,impossible,K,phi,pit,PMarkov,...
    Pout,rho,tau,tctp,theta0,tolInner,Vbsave)
% Approximate gradient wrt s


temp=Vbsave;
Vbsavetp=cell(1,tau);
for idx2 = 1:tau 
    Vbsavetp{idx2}=temp(:,:,:,idx2);
end
Nutp = size(tctp,1);
Pout0 = Pout;
step = 1e-4;
theta0(end) = theta0(end) + step;
s = theta0(end);

[Vb,~] = dreumInner5(H,impossible,K,tau,phi,PMarkov,rho,tctp,theta0,...
    tolInner,Vbsavetp);
[~,~,Pout,~] = dreumStationary_6(alpha,H,impossible,K,phi,pit,s,tau,Vb);

dPds = cell(1,tau);
for type = 1:tau
    for idx2 = 1:Nutp
        P0 = Pout0{type}(:,:,idx2);
        P = Pout{type}(:,:,idx2);
        dPds{type}(:,:,idx2) = (P - P0)./step;
    end
end
    